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Abstract 

A preconditioning scheme has been implemented into 
a three-dimensional viscous computational fluid dynamics 
code for turbomachine blade rows. The preconditioning 
allows the code, originally developed for simulating com- 
pressible flow fields, to be applied to nearly-incompress- 
ible, low Mach number flows. A brief description is given 
of the compressible Navier-Stokes equations for a rotating 
coordinate system, along with the preconditioning method 
employed. Details about the conservative formulation of 
artificial dissipation are provided, and different artificial 
dissipation schemes are discussed and compared. The 
preconditioned code was applied to a well-documented 
case involving the NASA large low-speed centrifugal 
compressor for which detailed experimental data are avail- 
able for comparison. Performance and flow field data are 
compared for the near-design operating point of the com- 
pressor, with generally good agreement between computa- 
tion and experiment. Further, significant differences 
between computational results for the different numerical 
implementations, revealing different levels of solution 
accuracy, are discussed. 


* Aerospace Engineer 

* Aerospace Engineer, Associate Fellow AIAA 

^Professor, Department of Mathematics, Senior Member AIAA 


Introduction 

The general engineering discipline of turbomachinery 
fluid dynamics involves a wide spectrum of practical 
devices and machines, many of which involve incompress- 
ible or nearly incompressible fluid flows. A list of such 
devices would include, among others, pumps, hydraulic 
turbines, propellers, low-speed fans and blowers, and low- 
speed experimental test rigs. Despite this abundance of 
incompressible-flow turbomachinery, and other incom- 
pressible flows in general, it is an interesting historical 
facet of computational fluid dynamics (CFD) that a great 
many computer codes for fluid dynamic simulation have 
been developed based on numerical algorithms intended 
only for compressible flows. One reason for this is that 
most state-of-the-art fluids engineering is being performed 
in the aeropropulsion industry where the fluid dynamics is 
mostly compressible, and where new technology is worth 
the high cost of expensive research and development, 
including CFD code development. 

It is well known that most compressible-flow CFD 
codes will not converge to an acceptable flow solution 
when the flow field Mach numbers become too low. Typi- 
cally this occurs somewhere below about Mach number 
0.1 where the gas becomes virtually incompressible; that 
is, the fluid density becomes nearly constant with changes 
in flow velocity. Within the past several years, however, 
methods for altering the compressible-flow numerical 
algorithm to allow convergence at very low Mach numbers 
have increasingly appeared in the literature [1-10]. These 
so-called preconditioning methods alter the eigenvalues of 
the system of compressible-flow equations so as to reduce, 
at low Mach numbers, the large disparity between the 
acoustic and convective wave speeds. Although it is possi- 
ble to avoid this problem altogether by developing the 
numerical algorithm to solve the incompressible flow 
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equations, the preconditioned compressible codes have 
several advantages and/or benefits, not the least of which 
is that many compressible -flow codes already exist for var- 
ious applications. It is relatively easy to add a precondi- 
tioning scheme to an already existing code. Other 
advantages include code versatility and the capability to 
directly simulate flow fields involving both compressible 
and incompressible flow regions. 

In this paper, a three-dimensional viscous CFD code 
originally developed for the simulation of the compress- 
ible flow field within a turbomachine blade row is 
considered. The code, designated RVC3D, solves a thin- 
layer formulation of the Reynolds-averaged three-dimen- 
sional Navier-Stokes equations [11, 12], and uses the 
Baldwin-Lomax algebraic turbulence model [13] to simu- 
late the effects of boundary layer turbulence. Precondi- 
tioning has been incorporated into the code by 
implementing the work of Turkei [3, 4], Radespiel and 
Turkel [7], Radespiel et al. [6], and Turkei et al. [10]. The 
paper describes the numerical method, including the pre- 
conditioning scheme, different artificial dissipation 
schemes, and some problems encountered. Computed 
results for a centrifugal compressor impeller are then pre- 
sented and compared to detailed experimental data. 

Governing Equations 

The governing equations which are numerically inte- 
grated using the RVC3D code are summarized below. 
Only the final thin-layer formulation of the viscous com- 
pressible equations as transformed to a generalized body- 
fitted coordinate system are presented here since a more 
comprehensive description can be found in the references 
for RVC3D [11, 12]. 

The Navier-Stokes equations are written for a Carte- 
sian (x, y, z) coordinate system rotating with angular veloc- 
ity O about the x-axis. The Cartesian equations have been 
mapped to a body-fitted (^, r|, Q coordinate system, sim- 
plified using the thin-layer approximation, and nondimen- 
sionalized by arbitrary reference quantities p 0 , a Q , and 

p 0 . Note that the ^-coordinate direction is roughly paral- 
lel to the blade surface and wraps around it, while the T|- 
direction is almost normal to the blade surface and the in- 
direction runs along the blade span. The resulting equa- 
tions can be summarized as follows: 

d t q + + d^G - Re~\d n F v + d^Gy) = fi 

( 1 ) 
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Details of the viscous flux vectors F v and G v are given 

in References 1 1 and 1 2. The absolute velocity compo- 
nents u, v and w point in the x , y, and z coordinate direc- 
tions, respectively, and the relative contravariant velocity 
components are given by 


V = 

v = ii*" + v' + 

W’ = + + 

where the relative velocity components are: 

u = u 
v' = v-&z 
w' = w + £ly 


(2) 


(3) 


then y ( = c p /c v ) is constant, and the energy and static 


Assuming an ideal gas with constant specific heats, 
1 Y ( = c p /c v ) is 
pressure are given by 

e = p 


1,2 2 2 / 
+ - ( U + V + w ) 


(4) 


Ly(Y-I) 2 

P = (Y-l)[e-ip(K 2 + v 2 + H> 2 ) (5) 

where the sonic velocity a is related to static pressure and 
density by the equation of state: 


2 y p 
a = — 


(6) 


The metric terms are defined by the following rela- 


tions: 


q = J q = J 
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n x C, 
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j y^-yfri 
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where the Jacobian J and its inverse J 1 are 


r ' s W! + W + Vft (8) 
•¥!V J i¥rW5 

The inlet (inflow) boundary condition for the above 
system of equations involves specification of the spanwise 
distributions of total pressure, total temperature, and abso- 
lute circumferential velocity component v 0 . The (non- 


preconditioned) up stream -running Riemann-invariant 

based on the meridional velocity component is extrapo- 
lated upstream from the interior to the inlet boundary. 

At the exit boundary four of the five conservation 
variables, namely p, pu, pv, and pw\ are extrapolated 
downstream to the boundary, and the static pressure is 
specified at the inner (hub) boundary and integrated in the 
spanwise direction using simple radial equilibrium: 


dp 

dr 


_ _ 2 

P^e 

r 


(9) 


where r is local radius and the overbar denotes a circum- 
ferentially-averaged quantity. The resulting local span- 
wise static pressures are then used either as constants 
(circumferentially-invariant exit static pressures), or as 
nominal values for calculating circumferential pressure 
distributions using a method described by Giles f 14]. 

Circumferentially -periodic boundaries between the 
blades are solved like interior points using a dummy grid 
line outside of the domain. At blade and endwall surfaces 
the no-slip, adiabatic wall boundary condition is used (for 
viscous simulations). 


Numerical Method 

The governing equations are discretized using a node- 
centered finite-difference scheme, with second-order cen- 
tral differences used throughout. The multistage Runge- 
Kutta scheme developed by Jameson, Schmidt, and Turkel 
[15] is used to advance the flow equations in time from an 
initial guess to a steady state. With the exception of the 5- 


stage scheme, the physical and artificial dissipation terms 
are calculated only at the first stage and then held constant 
for subsequent stages. Further details of the Runge-Kutta 
solution scheme can be found in References 1 1 and 12. 

The time step limit resulting from a three-dimensional 
linear stability analysis can be expressed as follows: 

* 

— — zr — =r — u (10) 

At ^ + A + Af^ +Af v 


* 

where A is the maximum Courant number for the partic- 
ular multistage scheme. The inverse one-dimensional 
time steps for each grid direction are equal to the largest 
local eigenvalue magnitudes (for the inviscid equations 
without preconditioning): 


where 
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(12) 


[T 2 2 


(13) 


The viscous time-step [16] contribution At 


-1 

v 


is given by 


At v = k ,^ a l + a l + a2 0 (,4) 

where the constant k f is assigned a value of 4.0. 

To accelerate convergence to a steady state, the maxi- 
mum permissible time step at each grid point is used, giv- 
ing a constant Courant number everywhere. To further 
accelerate convergence, implicit residual smoothing [17] 
is used to allow a time step typically two to three times 
larger than the stability limit given by Equation 10. Eigen- 
value scaling [18] of the residual smoothing coefficients, 
and of the artificial dissipation where applicable, is 
employed and greatly increases the accuracy and robust- 
ness of the numerical scheme. A blending function nearly 
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the same as one proposed by Kunz and Lakshminarayana 
[19] is used to blend the eigenvalues from all three direc- 
tions. 


Preconditioning 

The implementation of preconditioning into the 
RVC3D code is based on the work of Turkel [3, 4], Rade- 
spiel and Turkel [7], Radespiel et al. [6], and Turkel et al. 
[10], where generalized preconditioning schemes have 
been developed and summarized. The most recent publi- 
cations by the above authors present a generalization of 
the preconditioners originally introduced by Turkel [3, 4], 
and by Choi and Merkle [5], and they also present some 
further developments since then. In contrast to “standard” 
numerical schemes using a pseudo-compressibility 
approach, the preconditioners they describe not only 
accelerate the convergence of the numerical solution to a 
steady state, but also change the steady-state solution 
because of the choice of artificial viscosity terms. As dis- 
cussed in Reference 2, the standard pseudo-compressibil- 
ity schemes do not converge to the solution of the 
incompressible equations as the Mach number approaches 
zero. With proper preconditioning, however, the numeri- 
cal scheme does behave appropriately at low Mach 
numbers. The recent publications also show that precon- 
ditioning can be combined with well-known convergence 
acceleration techniques such as residual smoothing and 
multigrid. 

Only one of several possible schemes has been incor- 
porated into the RVC3D code, which is a scheme closely 
resembling the preconditioner of Choi and Merkel [5]. 
The scheme as implemented has been formulated to main- 
tain a conservative form of artificial dissipation. 

For convenience, Equation 1 is rewritten in a terse 
finite-difference form: 

A q = -A t[R,-(R v + D)] (15) 


Adi = -TAt 




( 17 ) 


(A^r 1 + A n r’v n + A ? r 1 v t > & 


where the solution vector q has been replaced by d) : 


d> — y (0 — 7 [p 9 /jj (18) 

and h is static enthalpy: 

h = a 2 /(y - 1) (19) 

The preconditioning matrix T and its inverse T _1 are 
given by 
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( 20 ) 
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-i 


i/p 2 

M/p 2 

v/P 2 

w/p 2 




0 0 0 0 
p 0 0 0 

0 p 0 0 

0 0 p 0 

p u pv pw p 


( 21 ) 


where Rj is the inviscid residual including source term, 
R v is the viscous residual, and D is the artificial dissipa- 
tion term. The artificial dissipation term can be expanded 
as follows: 

D = (A^ + A n V n + A ; V c )$ (16) 

where each coordinate-direction appears explicitly, and 
V^, , and V^, are operators defined later which 

depend on the type of artificial dissipation scheme. 
Denoting the preconditioning matrix by T, the precondi- 
tioned form of Equation 15 can be expressed as follows: 


2 2 2 2 

where q - u + v + w , and for the rotating coordinate 
system 

p 2 = min[max(?' 2 ,k q' 2 ,), a] 

P rej ( 22 ) 

,22,2,2 

q — u + v 

The parameters k p and q' re j are specified constants, with 
k p typically between 0.1 and 0.3 for viscous simulations, 
and q' re f being a nominal reference velocity in the rela- 
tive frame-of-reference. 

The largest eigenvalues for the resulting precondi- 
tioned (inviscid) equations are the following for each coor- 
dinate direction: 
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where 


'i + pT 

2 

v ay 


(24) 


and the relative contravariant velocity magnitudes \U'\ , 

| V'\ , and | W '\ , are obtained from Equation 12. 

The above eigenvalues are used to determine the time 
step and to directionally scale the implicit residual 
smoothing and the artificial dissipation in the same way as 
for the unpreconditioned scheme. 

Theoretical considerations suggest that modifications 
to the boundary condition formulations might be needed in 
conjunction with the preconditioning scheme [6, 10]. 
Good numerical convergence behavior was observed for 
all cases investigated, however, without such 

modifications. Therefore no changes were made to any of 
the boundary condition formulations. 


Artificial Dissipation Models 

Two basic artificial dissipation models are examined 
for use with the preconditioning scheme at low Mach 
numbers. The first model is that introduced by Jameson, 
Schmidt, and Turkel [15], and will be referred to as the 
baseline scheme. This scheme uses blended first- and 
third-order dissipation terms and is the one originally used 
in the RVC3D code as documented in the references [11, 
12]. The second model involves scaling the artificial dissi- 
pation by the local convective eigenvalue. At low Mach 
numbers (below Mach 0.5) this scheme is identical to the 
Convective Upwind Split Pressure (CUSP) scheme intro- 
duced fairly recently by Jameson [20], and subsequently 
modified by Tatsumi, Martinelli, and Jameson [21, 22]. 
The Symmetric Limited Positive (SLIP) flux limiter intro- 
duced by Jameson [20] is used in connection with this 
scheme, which for brevity will be referred to in this paper 
simply as the SLIP scheme. A variation of the SLIP 
scheme is also introduced, the ASymmetric Limited Posi- 
tive (ASLIP) scheme. For the low Mach number levels 
considered, the ASLIP scheme (convective-eigenvalue 
scaling with ASLIP flux limiter) appears to work slightly 


better than the SLIP scheme (convective-eigenvalue scal- 
ing with SLIP flux limiter). 

The baseline and the SLIP/ASLIP schemes can be 
described using the dissipation operators A^V^, A^V^, 

and in Equation 16 as the starting point. Since each 

directional operator is similar to the others, only the In- 
direction need be considered. Also, note that q every- 
where below is replaced by co when preconditioning is 
applied. 


Baseline Scheme 

In the baseline scheme, the dissipation term A^V^§ is 
finite-differenced according to 


A^V^ - ( J ^)/+i/2^; + 


1 / 2 ' 


.-1 


(25) 


(J CjO-. 


1/2 1 - 1/2 


where the /-indexing corresponds to the ^-direction and 
C ^ is a coefficient obtained from blended eigenvalues 
[19]: 


c % - h, 


1 + 


x„ + Xi 


% 


The dissipative flux vector 8 ;+ 1/2 is given by 


(26) 


^i + l/2 ~ V 2^i+l/2 

^4(^9, + 3/2 _ 2 A< 7 f+ 1/2 + _ i/ 2 ) 


(27) 


where A</ ( + 1/2 = q j+i -q i and the coefficients V 2 and 
V’ 4 are the following: 


V 2 = ji 2 mm(v f _,,v <1 v (+1 ,v ( + 2 ) 
V 4 = max(0,n 4 - V 2 ) 


(28) 


with 


= lP;+i— 2 P, + P/-il 

‘ xmn(p n ,\p i+l +2p i + p i _ l \) 


(29) 


The constant [i 2 scales a first-order artificial viscosity that 

is switched on at shocks detected by Equation 29. The 
denominator in Equation 29 is normally constant at the 
pressure p n , a fixed low pressure at either the grid inlet or 

exit, making the operator roughly symmetric across 
shocks. The more commonly applied term 
p i + t + 2 p i + /?,_]! is included to switch on the second- 
difference dissipation when the pressure becomes very 
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small, usually due to numerical problems. The constant _j_ 

ji 4 scales a uniform third-order artificial viscosity that is ? = max 0’ H — 2|A/ 


switched off at shocks by V 2 • Typically, values for p 2 
and p 4 are 1/4 and 1/32 , respectively. 

SLIP/ASLIP Scheme 


S = min(l, 


S + C c 


( 34 ) 


In the SLIP/ASLIP scheme, the dissipation term 
is finite-differenced according to 


= {7 ~ 


1/2 ‘ 


,-1 


(30) 


<' °5><- 


l v i- 1/2 


where is given by Equation 13, and the dissipative flux 
vector 5- + 1/2 is given by 


*>1 + 1/2 * 2 aa/+1/2 ^* ^ 

The subscripts R and L denote right and left states, respec- 
tively, for the state vector q. The left and right states are 
determined by either the SLIP or the ASLIP flux limiter, 
as described later. The coefficient a is given by 


|Af 


'd 


a = q 


M 


>2\ 


q eo + 


\M'^ > q eo 
jM'^l < q eo 


(32) 


q eo J 


where the relative contravariant Mach number is 




V 


V’r + V'l 


1/2 2 ^) <+ l/2 


(33) 


When preconditioning is applied, q R - q L is replaced by 
r.| 1 / 2 ( 0 )^ - co L ) in Equation 31. 

The parameters and q rescale the artificial dissi- 
pation near the wall for viscous flows, for reasons dis- 
cussed in Reference 23. For convenience and 
computational efficiency, the coefficients obtained from 
blended eigenvalues (see Equation 26) are used: 


The parameter q as given above is used only in the In- 
direction, with q+ being set identically to one for the T|- 
and indirections (where physical viscosity terms exist). 
The parameters e Q is defined by 

E 0 = K M \ef ( 35 ) 

where k v and M' re f are specified constants, with k y typi- 
cally given a value of 0.05 or less for viscous simulations, 
and M' re f being a nominal reference Mach number in the 

relative frame-of-reference. 

As alluded to earlier, either the SLIP or ASLIP flux 
limiter is used to determine the left and right states. This 
is done by extrapolation from nearby data, subject to a 
limiter to preserve monotonicity. This approach is similar 
to van Leer’s MUSCL scheme [24] and provides first 
order dissipation near “discontinuities”, that is, non- 
smooth regions, and third-order dissipation in smooth 
regions of the flow field. 

From the work of Swanson, Radespiel, and Turkel 
[23], and considering the right and left states in the dissi- 
pative flux vector + 1/2 , a slightly generalized form of 
the SLIP limiter can be written as follows: 

U = ^/ + 2 L ^i + 3/2’ A ^+l/2’ A< 7(-l/2) 

x (36) 

Qr = 4| + 1 “2^ A ^ , + 3/2 ’ ^i+l/2’ A ^/-l/ 2 ) 

where A^ J + 1/2 = <7/ +]-“#; , and the operator L is 
defined as 


A q, Aq j = 

(1 

with the switch S as 
S^A<? + , Aq j = 


(a/, Aq j 

+ I 2 


A <5 

+ 

-Aq 

X 

a/ 

+ 

1 

< 

1 

1 

I + £ 


(37) 


( 38 ) 
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The exponent % is a positive number given a value of 2.0 in 
the RVC3D code. The constant fi 4 scales the third-order 
dissipation such that in smooth regions of the flow 
(5-1): 


H 4a 3 

~ 2 A + 1/2 


(39) 


Typically |i 4 is specified to have a value of 1 .0. Note that 

whenever A q and A q have the opposite sign then 

5 — 0 , and where the flow is smooth then 5 — 1 . The 
parameter £ in Equation 38 is a threshold to prevent the 

denominator from getting close to zero when both A q 

and A q are small. 

The ASLIP limiter is similar to the SLIP limiter 
except that the operator L is split into two parts: 


L^Aq+ y Aq, Aq (a? ,A q^j+^L (a<? , A q 

(40) 

where 

L + ^A^ + , Aq'j = S^Aq + , A<?j 

f (l - 2|i 4 )A q + H 4 (a<? + + Aq ) 

(41) 

L , A<?^ = S^Aq , Aq j 

j/ 1 - 2|i 4 )A q + H 4 (a<? + Ag) 


The switch 5 is defined the same as above. Note that in 
smooth regions of the flow the ASLIP limiter is virtually 
the same as the SLIP limiter, but in non-smooth regions 
the ASLIP limiter becomes asymmetric about i + 1/2 due 

to the independent action of the two halves, L and L , 
of operator L. 


Comparison of Baseline and SLIP/ASLIP Schemes 

A rigorous comparison of the above schemes is well 
beyond the scope of this paper, but a couple of preliminary 
comparisons are provided to demonstrate significant dif- 
ferences between the schemes. In particular, results for 
the nearly-incompressible laminar boundary- layer flow 
over a flat plate are presented for two conditions, 
freestream Mach numbers M ^ - 0.30 and M ^ — 0.05 , 

respectively. For both flows the same grid was used and 
the Reynolds number based on plate length was 

Re L = l.OxlO 6 . Note that instead of RVC3D, a two- 

dimensional CFD code [25] was used for these simula- 
tions, but the numerical formulations are very similar for 
both codes. Regarding the grid (not shown due to space 
limitations; 41 nodes normal to wall; 145 nodes along 
plate), it should be noted that grid characteristics such as 
coarseness, wall spacing, and wall stretching, were 
selected to be representative of typical RVC3D computa- 
tions, and that these characteristics substantially influence 
the numerical solutions. Typical values for the numerical 
parameters were used as well (p 4 = 1/32 for the base- 
line scheme, and |X 4 = 1 for the SLIP and ASLIP 

schemes). Finally, it should be emphasized that the com- 
parisons made here apply only to Mach numbers less than 
0.5. 

Although the flow over the entire plate length was 
simulated in each case, it is sufficient to compare the 

results at only one plate location, Re x = 5.0xl0 5 . Simi- 
lar differences were observed at other plate locations. 

The results for freestream Mach number M ^ = 0.30 
are shown in Figure 1 where the normalized errors 
Au/u ^ in the w-velocity profiles are compared for the 

baseline scheme (with and without preconditioning) and 
the SLIP and ASLIP schemes (only with 
preconditioning). A local error of zero corresponds to a 
local u-velocity which is equal to that of the “exact” Bla- 
sius solution. The dimensionless distance from the wall is 

the standard similarity parameter r\ = (y/x)jRe x /2 . 

The comparisons in Figure 1 reveal a relatively large 
difference in accuracy between the baseline and the SLIP/ 
ASLIP schemes, with the latter performing considerably 
better. It is noteworthy that the baseline scheme is 
improved substantially by the use of preconditioning, 
which in this case reduces the artificial dissipation (see 
Equations 25 and 26) due the to rescaling of eigenvalues 
(compare Equations 1 1 and 23). Although outside the 
scale of the graph, the maximum error of around -7.5 per- 
cent for the baseline scheme without preconditioning is 
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Normalized Velocity Error, Au/u^, Percent 

Figure 1 — Error in computed solutions of laminar Blaslus 
velocity profiles using the baseline, SLIP, and 
ASLIP schemes (M^ = 0.30) 


reduced down to about -3 percent with preconditioning. 
This result illustrates the fact that at low Mach numbers 
the preconditioning methods being implemented favorably 
alter the steady-state solution, as noted earlier. The SLIP 
and ASLIP schemes both perform well, the ASLIP scheme 
being slightly more accurate than the SLIP. This small 
difference has been confirmed in general for a wide range 
of (undocumented) practical low Mach number simula- 
tions for various turbomachinery cases. It should be 
noted, however, that the ASLIP flux limiter has been 
observed to work well only for low Mach number flows, 
whereas the SLIP flux limiter has a far wider range of 
application [20-22]. Finally, it can be seen in Figure 1 that 
the SLIP and ASLIP schemes completely eliminate the 
velocity “overshoot” oscillation at the outer edge of the 
layer, an effect commonly observed with respect to the 
baseline scheme at virtually all Mach number conditions. 
For a transonic airfoil example involving such an over- 
shoot and its elimination by using the CUSP scheme with 
the SLIP flux limiter [20-22], see Reference 26. 

Laminar boundary layer results similar to those above 
can be seen in Figure 2 for the freestream Mach number 
M ^ = 0.05 . In this figure, however, no curve is included 

for the baseline scheme without preconditioning since the 
low freestream Mach number makes that computation 
impractical. Qualitatively the results are the same as for 



Normalized Velocity Error, Au/u^ , Percent 

Figure 2 — Error in computed solutions of laminar Blaslus 
velocity profiles using the baseline, SLIP, and 
ASLIP schemes (M m = 0.05) 

M ^ = 0.30 , although quantitatively the near-wall errors 

for the SLIP and ASLIP schemes are somewhat larger. 

The above comparisons suggest that improved accu- 
racy might be obtained in general by using the SLIP or 
ASLIP scheme rather than the baseline scheme. However, 
there are some disadvantages worth mentioning. A pri- 
mary disadvantage is the direct increased computational 
overhead. For the 4-stage Runge-Kutta scheme with 
implicit residual smoothing, the overall CPU-time per iter- 
ation increases by about 26 percent for a typical low -speed 
viscous simulation, while for the 2-stage Runge-Kutta 
scheme the increase is about 50 percent. Note that, as 
implemented in RVC3D, each scheme applies the artificial 
dissipation only once per time-step. Another notable dis- 
advantage is that the SLIP and ASLIP schemes sometimes 
reduce the code stability with 4-stage Runge-Kutta time- 
stepping, although the most frequent adverse effect is just 
poor convergence behavior. The exact cause of this effect 
has not been determined, but it is thought that a mu hi grid 
scheme (not in RVC3D) might eliminate or reduce it while 
also improving code stability. Interestingly, good stability 
and convergence behavior is always achieved when using 
the SLIP and ASLIP schemes in connection with 2-stage 
Runge-Kutta time-stepping where, with implicit residual 
smoothing, CFL numbers of 2.5 to 3.0 can be routinely 
used. 
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Numerical Source Terms 

Application of the RVC3D code to low Mach number 
flows revealed that some approximations related to the 
finite -differencing were inadequate, even though no signif- 
icant effects were encountered at higher Mach numbers 
where the code has been routinely applied. As a result, 
two basic modifications were necessary in order to avoid 
substantial numerical sources and get good convergence 
behavior with undistorted flow field solutions. First, the 
expressions for the finite-difference metrics had to be 
replaced with more rigorously formulated (and conse- 
quently more complicated) ones. Second, a corrective 
source term had to be added to account for errors intro- 
duced by coordinate-system rotation. 

Finite-Difference Metrics 

Until recently the RVC3D code used second-order 
central differences to evaluate the metric terms in 
Equation 7. However, as was pointed out years ago by 
Pulliam et al. [27], in three dimensions this scheme can 
introduce a freestream error, although it is noteworthy that 
Reference 27 mentions no discernible effects from ignor- 
ing the numerical source error other than the inability to 
maintain the free stream. At low enough Mach numbers, 
however, the resulting numerical sources are large enough 
compared to the physical flux terms to produce excessive 
flow field distortions, as well as hinder or prevent numeri- 
cal convergence. 

Metric finite-difference formulas based on cell-face 
projected areas have been in existence for some time, and 
can be found clearly presented in Reference 28. For con- 
venience, the numerical formulas are also given here: 

2 J = A y(0, 1, l) A z(0,-l, l) _A y(0,-l, l) A z(0, 1, 1) 

2 J = A z(0, 1, 1) A *(0,-1, 1)~ A z(0,-1, 1) A *(0, 1, 1) 

2 J = A *(0, 1, 1 ) A v(0, -1,1)" A *(0, -1, l) A y(0, 1, 1 ) 

2 J *1* = A y(l, 0, 1 ) A z(l, 0, -1) — A y(l, 0, -l) A z( 1, 0, 1 ) 

2 J % = A z(1,0,1) A x(1,0,-1)“ A z(1,0,-1) A jc(1,0, 1) 

2 J *1* = A *(l,0, 1 ) A >(1, 0, -1) “ A x(l,0,-l) A y(l,a 1) 

2 J C* = A y(l, l,0) A z(-l, 1,0) “ A y(-1, l,0) A z(l, 1,0) 

2 J £>y = A z (l, l,0) A x(-l, l,0)“ A z(-l, 1,0 ) A jc(1, 1,0) 

2 J = A jc(l f 1, 0) A v(-l, 1, 0) “ A jr(-1, 1, 0) A v(l, 1, 0) 

(42) 


where the difference operator is defined as 

A x(I,J,K ) = 2 ( X i + IJ + J, k + K~ Xi -IJ- J , k -K^ 

Note that the indices y, and k correspond to the coordi- 
nate directions T|, and £, respectively, and that these 
expressions for the metrics are fully conservative. 

The metric Jacobian J is calculated as before, using 
central-finite-differences, since it does not significantly 
influence the steady-state solution. However, unlike the 
case of a stationary coordinate system there is some (negli- 
gible) influence due to the rotation source term on the 
right-hand-side of Equation 1 . 

Coordinate-System Rotation 

The numerical source error due to coordinate-system 
rotation can be examined analytically by considering a 
special case; namely, that of steady uniform flow in the x- 
direction, in which case the velocity component u is a con- 
stant, the components v and w are zero, and all flow deriv- 
atives are zero. The continuity, x-momentum, and energy 
equations (see Equation 1) can then be reduced to the fol- 
lowing expression: 

0 = n [d^-r'f, y z + J~ X t, z y) + 

a n (-.r 1 v + - / ” 1 Tvv)+ < 43 > 

The right-hand side of Equation 43 is mathematically 
equivalent to 

A By dz) 

which is identically zero for any non-zero angular velocity 
£1. For the numerical formulation the right-hand side is 
not exactly zero, however, which produces a freestream 
error. 

Using the result of Equation 43, and denoting the 
right-hand side as a simple and straightforward 

source term correction can be applied in Equation 1. In 
particular, an additional source term k can be included to 
exactly cancel the freestream error: 

d,q + + 3 ; G - Re~' 0/ v + ^G v ) ^ 

= // + k 

where 

k = S^q (45) 


NASA TM-1 13120 


9 



Note that for efficiency, k is only evaluated at the first 
stage of the Runge-Kutta scheme and then held constant 
along with the artificial dissipation and the viscous resid- 
ual. 


Numerical Results 

The preconditioned RVC3D code has been applied to 
a well-documented experimental test case involving a 
large low-speed centrifugal compressor (LSCC) impeller 
tested at the NASA Lewis Research Center. After describ- 
ing the experimental test facility along with the CFD grid, 
several computational details are discussed, followed by 
comparisons of computational and experimental results. 

Low-Speed Centrifugal Compressor Impeller 

The LSCC is a large-scale low-speed compressor test 
rig designed to duplicate the essential flow field and flow 
physics of a high-speed centrifugal compressor impeller 
[29-31]. Due to its large size and low speed, extensive 
aerodynamic instrumentation of the impeller was made 
feasible, as were detailed flow field measurements using a 
laser anemometer system [31]. 

In order to provide an illustration of the test impeller, 
a coarse rendering of the CFD grid on the impeller hub 
and blade surfaces is shown in Figure 3. The impeller has 
20 full blades with 55 degrees backsweep and a design tip 
speed of 153 m/s (1917 rpm). The inlet diameter is 
0.870 m with an inlet blade height of 0.21 8 m, and the exit 
diameter is 1.524 m with an exit blade height of 0.141 m. 
The tip clearance between shroud and impeller is uniform 
at 2.54 mm, which is 1.8 percent of blade height at the 
impeller exit. Fillet radii at the impeller hub are uniform 



at 9.525 mm. A vaneless radial diffuser is located down- 
stream of the impeller. 

Although Figure 3 shows all 20 impeller blades, the 
computational grid for RVC3D has only one blade (and 
blade channel), with blade-to-blade flow periodicity 
enforced by the boundary conditions. The RVC3D code 
uses C-type grids [1 1,12], and the same grid was used for 
all computational results presented below. The blade-to- 
blade grid near midspan is shown in Figure 4 which 
includes expanded views for the blade leading- and trail- 
ing-edges. The dimension of the grid in the ^-direction 
(wrapping around the blade) is 225, while in the i\ -direc- 
tion (normal to the blade surface) the grid has a dimension 
of 33, giving a total of 65 nodes blade-to-blade. In the in- 
direction (spanwise, hub-to-shroud), the grid dimension is 
65 with more clustering near the shroud for better resolu- 
tion of the blade tip clearance flow. As shown in Figure 5, 
the grid was stretched into the clearance gap region over 
the blade tip so as to provide for a direct simulation of the 
tip clearance flow, but without having to use a multi-block 
flow solver. The blade-hub fillets were not included in the 



Figure 3 — LSCC impeller with coarse rendering of compu- 
tational grid on hub and blade surfaces 


Figure 4 — Near midspan blade-to-blade view of computa- 
tional grid for LSCC impeller 
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Figure 5 — Computational grid stretching into shroud 
clearance gap over blade tip 


computational grid geometry since they are thought to 
have only a small influence on the impeller performance 
and flow field. 

Near-Design Operating Condition 

The compressor operating condition of interest for the 
present study is a near-design point for which extensive 
measurement data were obtained. The corrected flow rate 
at the near-design point was 30.0 kg/s with the impeller 
rotating at 1862 rpm. Based on aerodynamic probe mea- 
surements, the overall impeller performance at that point 
has been reported [31] to be a total -pres sure ratio of 1.141 
with an adiabatic efficiency of 0.922. 

Inlet Duct Flow Simulation 

To more accurately simulate the flow entering the 
impeller, while also eliminating several approximations 
often made at a rotor grid inlet boundary, an axi symmetric 
flow simulation for the entire inlet duct upstream of the 
impeller was performed by “loosely” coupling the duct 
fluid dynamic simulation to the impeller simulation. The 
coupling was implemented through the appropriate bound- 
ary conditions, and at the level of the computer operating 
system. Specifically, for every 150 iterations of the impel- 
ler simulation, a converged axisymmetric duct simulation 
was performed using the circumferentially-averaged span- 
wise static pressure distribution from the impeller grid 
inlet. The viscous two-dimensional CFD code used for 


simulating the axisymmetric duct flow is described in 
Reference 25. The resulting spanwise distributions of 
total pressure, total temperature, circumferential velocity, 
and meridional flow angle at the duct exit were then 
passed to the impeller grid inlet for the next block of 150 
iterations. Good overall convergence between the duct 
and impeller solutions were obtained by under-relaxing 
the duct exit static pressure update at the beginning of each 
duct simulation. 

Resulting velocity contours for the duct and impeller 
solutions as computed using preconditioning and the 
ASLEP artificial dissipation scheme (both solutions) are 
shown in Figure 6, where the velocities in the LSCC 
impeller are circumferential ly-averaged and shown in the 
absolute frame-of-reference. Although mass-flow conti- 
nuity was not strictly enforced between the duct and the 
impeller, the integrated values are nearly the same at 29.96 
and 29.94 kg/s, respectively, which are within 0.2 percent 
of the measured flow rate. 

General Discussion of Impeller Simulations 

Several numerical simulations were performed for the 
LSCC impeller. It will be helpful to delineate them before 
proceeding to detailed comparisons of computational and 
experimental results. Numerical convergence issues will 
also be addressed in this context. 

Two separate but related numerical scheme imple- 
mentations are being examined in this study: namely, pre- 
conditioning and artificial dissipation. Use or non-use of 
the preconditioning with any of the artificial dissipation 
schemes discussed earlier leads to many possible combi- 
nations which could be compared. Due to scope and space 
limitations, however, only three are considered in this 
study: 

• The baseline artificial dissipation scheme with no 
preconditioning. This case is possible due to suffi- 
cient compressibility of the flow being simulated, and 
allows comparisons of results computed with and 
without the preconditioner. 

• The baseline artificial dissipation scheme with pre- 
conditioning. 

• The ASLEP artificial dissipation scheme with 
preconditioning. The ASLEP was selected over the 
SLIP scheme since it is thought to be the most accu- 
rate at the low Mach number levels under 
consideration. However, the SLIP scheme generally 
produces results similar to those of the ASLEP 
scheme. 

For the baseline case with preconditioning, as for the 
ASLIP case, the solution was obtained with coupling to 
the inlet duct flow. However, the duct solution is virtually 
the same as that shown in Figure 6. The baseline case 
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Figure 6 — Inlet duct and circumferentially-averaged LSCC impeller flow fields 


with no preconditioning was solved without inlet cou- 
pling, using instead the inlet duct solution from the pre- 
conditioned baseline case. 

In general, preconditioning can be expected to accel- 
erate solution convergence for low Mach number flows, 
and if the Mach levels are low enough should allow a con- 
verged solution where none is possible without it. For the 
present case a convergence-history comparison is possible 
since the simulation converges without preconditioning. 
For the baseline artificial dissipation scheme (p 4 = 1/32 

and \i 2 = 1 /4 ), the benefit provided by preconditioning 

is shown in Figures 7 and 8. The normalized RMS-resid- 
ual history in Figure 7 reveals a large acceleration in con- 
vergence over the first 2000 iterations. The accompanying 
variations in global mass-flow error and overall total-pres- 
sure ratio in Figure 8 are perhaps more revealing in that 
they show a dramatic acceleration from the precondition- 
ing, such that by 1000 iterations the large oscillations in 
global mass-flow error and overall total-pressure ratio 
have been greatly diminished. It should be noted that both 
of these simulations were performed without coupling to 
the inlet duct in order to maintain a simple, direct 
comparison. They were executed using the 4-stage 
Runga-Kutta time-integration with implicit residual 
smoothing and a CFL number of 5.6. The computational 
overhead incurred by using the preconditioner for this sim- 
ulation was determined on the CRAY C-90 to be about 22 
percent, but as can be seen in the figures the convergence 
rate has been more than doubled as a result. Actual execu- 
tion speeds for the non-preconditioned and the precondi- 
tioned solutions on the CRAY C-90 (using a single CPU) 
were 2.17 and 2.64 sec/iteration, respectively. 

As noted earlier, the use of the ASLIP (or SLIP) artifi- 
cial dissipation scheme with the 4-stage Runga-Kutta 
time-integration often results in poor convergence 
behavior. Although the effect seems more associated with 



Number of Iterations 


Figure 7 — Comparison of residual histories for simula- 
tions with and without preconditioning (base- 
line artificial dissipation scheme) 


the residual history rather than the overall flow solution — 
that is, a relatively small local conveigence problem 
causes the residuals to “hang up” — the problem is com- 
pletely avoided by using a 2-stage Runga-Kutta scheme 
where, with implicit residual smoothing, a CFL number of 
3.0 can typically be used. In general, the use of the 2- 
stage Runge-Kutta solver requires roughly 50 percent 
more CPU time (by requiring more iterations) to achieve a 
converged solution. For the LSCC impeller on the CRAY 
C-90, the execution speed with the ASLIP scheme and 
preconditioning was 2.67 sec/iteration for the 2-stage 


NASA TM-1 13120 


12 





Figure 6 — Comparison of convergence histories for simu- 
lations with and without preconditioning (base- 
line artificial dissipation scheme) 


solver, and about 3.31 sec/iteration with the 4-stage 
solver. Note that all of the ASLEP results presented below 
were obtained using the 2-stage solver, whereas the base- 
line cases were run using the 4-stage solver. 

Detailed Comparisons of Results 

The computational results are compared with experi- 
mental data and with each other at several locations within 
or near the LSCC impeller. The experimental measure- 
ment locations of interest are shown in Figure 9 where the 
measurement station J-numbers are indicated for later 
reference. Computational results are first compared with 
experimental data at the two measurement stations 
upstream of the impeller (7 = 23 and 48) in order to estab- 
lish that the computed and measured impeller inflow con- 
ditions were in fact nearly the same. 


Spanwise distributions of meridional velocity V m 

nondimensionalized by impeller tip speed U tip , and 

meridional flow angle a. m are compared at station 23 in 

Figure 10, and at station 48 in Figure 11. As Figure 9 
reveals, the computed solution at station 23 is within the 
inlet-duct domain since the grid interface (between the 
duct grid and the impeller grid) is downstream of 
station 23. Conversely, station 48 is located within the 
impeller-grid domain and is just upstream of the impeller 
leading edge. Since all three computed solutions are simi- 
lar upstream of the impeller, only one computed result 
(ASLIP scheme) is shown for each station. The experi- 
mental data were acquired with a laser Doppler 
velocimeter (LDV), and circumferentially -distributed 
measurements have been circumferentially mass-averaged 
to obtain the results graphed. Generally close agreement 
between the numerical and experimental results can be 
seen in the comparisons, confirming the accuracy of the 
computed impeller inflow. 

As the flow approaches and enters the impeller it 
accelerates and therefore undergoes some drop in static 
pressure. Within the impeller, however, there is a steady 
rise in both static and total pressure as the flow moves 


Rotor 

Trailing 

Edge 



Figure 9 — Experimental measurement station locations 
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Figure 10 — Comparison of computed and measured span- 
wise profiles at Impeller inlet station 23 



Figure 1 1 — Comparison of computed and measured span- 
wise profiles at Impeller inlet station 48 


through it and the impeller does work on the fluid. These 
features are readily apparent in Figure 1 2 which compares 
the measured steady-state shroud (tip) static pressure dis- 
tribution through the impeller, to the computed and cir- 
cumferentially-averaged shroud static pressure 
distribution. The pressures have been normalized by the 
inlet total pressure p 0 , and the results are graphed as a 

function of meridional distance m (from the impeller lead- 
ing edge) normalized by the total meridional distance 
through the impeller blades . As with the inlet 

profiles, the differences between the various computed 
solutions is negligibly small so only one of them need be 
presented. It can be seen that the agreement between com- 
putation and experiment is very good, except perhaps at 
the impeller exit (m/ Am^ ^ 100 percent) where some 
differences are apparent. 

Considerably more information about the pressure 
field within the impeller can be obtained from the static 
pressure distributions for the impeller blade surfaces. 
Blade surface pressure distributions for six different span- 
wise locations, namely 5, 20, 49, 79, 93, and 97 percent 
span from the hub, are compared in Figure 13. Again, the 
differences between the three computed solutions are 
small enough to be ignored in this figure. As can be seen, 
the agreement between computation and experiment is 
very good, except on the pressure surface (PS) near the 
blade trailing edge ( m / ^ > 90 percent) where sig- 



Percent Meridional Distance, m /Am kl _ 

blade 

Figure 12 — Comparison of computed and measured 
shroud static pressure distributions 
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Blade Surface Static Pressure, p/p 



Percent Meridional Distance from Leading Edge, m/Am bfgde 

Figure 13 — Comparison of computed and measured blade static pressure distributions at several spanwise locations 
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nificant differences exist for most span wise locations. 
Notice that at each spanwise location both the suction sur- 
face (SS) and pressure surface (PS) distributions are 
graphed together, which conveys useful information about 
the blade static pressure loading, the fundamental means 
by which the impeller does work on the fluid. The greater 
the area between the SS and PS curves, the greater the 
blade loading. With this in mind, examination of the 
graphs reveals that the experimental impeller has substan- 
tially more loading near the trailing edge than indicated by 
the numerical simulations. It can therefore be anticipated 
that the overall work input and total -pressure rise obtained 
from the computed results will be less than measured. 
Later it will be shown that this is actually the case. 

In terms of the static pressure and inlet profile results 
presented so far, only results from one of the computed 
solutions (the ASLIP solution) were provided since, as 
noted, only minor differences existed between the three 
computed results. This does not imply, however, that only 
minor differences exist in general between the three 
solutions. In fact, significant differences do exist as will 
become apparent in much of what follows. 

Contour plots of the blade-to-blade flow field at the 49 
percent spanwise location from the hub are shown in 
Figure 14 for the three different computed cases (baseline 
with no preconditioning, baseline with preconditioning, 
ASLIP with preconditioning). Two sets of contour plots 
are compared: relative velocity V'/« 0 on the left, and rel- 
ative total pressure p f '//? 0 on the right. Note that the 

contour levels are not labeled in this figure since only a 
qualitative comparison is to be made. Close examination 
reveals that there are clear differences between the three 
solutions and, furthermore, that the differences are of such 
a character that a reasonable inference with regard to rela- 
tive accuracy can be made. 

Inspection of the relative total pressure contours in 
Figure 14 while taking into consideration certain basic 
fluid dynamic principles leads to a few useful 
conclusions. First, it can be pointed out that apart from 
any viscous flow effects the blade-to-blade contours should 
be straight, vertical lines in the figure. This follows from 
the fact that for a homentropic inlet flow without swirl the 
relative total enthalpy, and thus the relative total pressure, 
should depend only on radius; that is, the “rothalpy” is 
constant along an inviscid streamline. Viscous regions of 
the flow field such as the blade and endwall boundary lay- 
ers, blade wakes, and any viscous secondary flows are 
expected to produce distortions in the contour shape. It 
can be seen in the figure that only the ASLIP solution 
nearly achieves the expected result — straight vertical 
lines between blade boundary layers — whereas the base- 
line scheme without preconditioning departs substantially 
from it. The corresponding differences in the midspan 


velocity field (Figure 14, left) are significant, especially in 
the exit region of the impeller. With preconditioning the 
baseline scheme is significantly better than without it, but 
noticeable errors in the total pressure contour shapes can 
still be seen. Also noticeable in the two baseline solutions 
but absent in the ASLIP solution is a total pressure oscilla- 
tion (overshoot) at the outer edge of the blade boundary 
layers. The velocity contours also exhibit this overshoot, 
although it is less apparent. Recall that this behavior was 
shown and discussed earlier for the laminar flat plate 
boundary layer solutions (see Figures 1 and 2). In fact, the 
general trend in accuracy for the three LSCC impeller 
solutions parallels exactly that presented earlier for the flat 
plate. 

Detailed comparisons of the computed and measured 
results downstream of the impeller at station 178 (see 
Figure 9) are presented next, followed by some compari- 
sons at stations within the impeller. At this juncture it 
might be pointed out that in the aft half of the impeller 
strong viscous secondary flows develop and that these 
flows reside mostly in the outer one-third (shroud side) of 
the span while within the impeller. Downstream of the 
impeller, spanwise mixing and migration begin to redis- 
tribute the flow. 

Spanwise distributions of steady- state total pressure 
P/Pq and static pressure p/p 0 at station 178 are com- 
pared in Figure 15. Note that the hub and shroud corre- 
spond to the right and left sides, respectively, of each 
graph. The experimental data in this figure were obtained 
using slow-response pneumatic probes, and the computa- 
tional results represent circumferentially mass-averaged 
solution data. Closely related to the total pressure distri- 

_ 2 

butions are the work-coefficient A (UV^)/U tip distribu- 
tions presented in Figure 16 where results for station 165 
(see Figure 9) are included below those for station 178. 
Again, the computational results represent circumferential 
mass-averages, as do the LDV experimental data. The 
experimental “pneumatic probe” results in Figure 16 were 
calculated from measured total pressures using the well- 
known Euler work equation, as described in Reference 31 . 

The static pressure distributions in Figure 15 (bottom) 
were determined primarily by the static pressure condi- 
tions imposed downstream, either by the throttle setting in 
the experiment or by the grid exit static pressure distribu- 
tion in the computations. Slightly different grid exit static 
pressures, constant over the span, were specified for the 
three CFD cases so as to obtain the same mass flow rate. 
As can be seen in the graph, the experimental data indicate 
a slight positive pressure gradient from shroud-to-hub — 
which is assumed real, but could be a measurement error 
due to probe immersion blockage effects — but the magni- 
tude is small enough to be neglected. 
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Figure 15 — Comparison of computed and measured span- 
wise total and static pressure distributions at 
station 178 


The total pressure distributions in Figure 15 (top) 
reveal significant differences between the computed 
results, all three of which differ quite substantially from 
the experimental results. The work coefficient distribu- 
tions in Figure 16 (top) show virtually the same trends and 
differences, which strongly implies that the differences in 
total pressure are due to differences in impeller work 
rather than in impeller total-pressure loss. As was pointed 
out earlier, the higher level of experimental impeller work 
is to be expected in view of the differences in blade (static 
pressure) loading (see Figure 13). Recall from Figure 13 
that the experimental blade loading over the last 5 to 10 
percent of blade chord (near the trailing edge) was discem- 
ibly higher than computed. In fact, closer examination of 
the pressure distributions in Figure 13 makes it clear that 
the simulated impeller undergoes a rapid aerodynamic 
unloading, even a negative loading, as the trailing edge is 
approached. This unloading is due mostly to an abrupt 
drop in static pressure on the blade pressure side. The 
experimental pressure surface data, however, do not show 
this. 




Percent Passage Height 

Figure 16 — Comparison of computed and measured span- 
wise work-coefficient distributions at stations 
178 and 165 


Further substantiation of the above discussion is pro- 
vided by the spanwise work coefficient distributions at 
station 165 (see Figure 16, bottom), which is located in the 
impeller at about 94 percent of blade chord. At that loca- 
tion it can be seen that the agreement between computa- 
tion and experiment is much better, except around 80 
percent of passage height (span) where the influence of a 
viscous, secondary-flow tip-clearance vortex is strong. 
More spanwise comparisons at stations 178 and 165 

are given in Figure 17; radial velocity V r /U fi on the top 

half, and relative flow angle p on the bottom half. Note 
again, that the quantities graphed are circumferential-aver- 
ages, and that the flow angles were calculated using aver- 
aged circumferential and radial velocities: 

P = tan -1 ( V e '/V r ) (46) 
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ASLiP 



Percent Passage Height 

Figure 17 — Comparison of computed and measured down* 
stream work-coefficient distributions 


It might be added that at stations 165 and 178, the radial 
velocity is approximately equal to the through flow veloc- 
ity since the flow path is nearly radial. 

As can be seen in Figure 17, the computed radial 
velocities at station 165 are in good agreement with each 
other and with the experimental data. At station 178 the 
agreement in radial velocity between computation and 
experiment is less favorable, but still good, and there are 
larger differences between the computed solutions. The 
relative flow angle distributions for station 165, like the 


radial and work coefficient (absolute tangential velocity) 
distributions, are also in reasonable agreement, except 
near 80 percent span (compare Figure 16, bottom). The 
differences between the computed and measured angle 
distributions at station 178 are merely a consequence of 
the differences in impeller work since the radial velocities 
are in fair agreement. 

An attempt to understand and explain the underlying 
reason(s) for the discrepancies between computation and 
experiment at station 178 is well beyond the intent and 
scope of this study. Nevertheless, some further description 
and explanation should provide a better perspective to the 
comparisons being made, and might also help to stimulate 
more insight for future investigations. 

An overall picture of what is happening at the impel- 
ler exit can be seen by considering a nominal velocity vec- 
tor triangle for the exit flow. This consideration need only 
be qualitative, as in Figure 1 8 (top) showing the basic dif- 
ferences between the experimentally and computationally 
derived triangles. The triangles in the figure have been 
drawn such that the absolute tangential velocity V 0 
(related to impeller work; see Figure 16, top) is larger for 
the experimental flow, and the radial velocity V r is about 
the same for computation and experiment (see Figure 17, 


4 



Figure 18— Qualitative comparison of exit velocity vector 
triangles and impeller trailing-edge flows 
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top). As a result, the relative tangential flow velocity V 0 ' 

and relative flow angle P are smaller for the experimental 
flow (see Figure 17, middle bottom). Since the differences 
illustrated in Figure 18 are generated in just the aft 5 to 10 
percent of the impeller, and since the computed pressure- 
surface pressures on the blade drop rather abruptly near 
the trai ling-edge whereas the experimental pressures do 
not (see Figure 13), it can be inferred that the computed 
flow experiences a trailing-edge deflection away from the 
direction of impeller rotation (^/-direction) while the 
experimental flow does not, as illustrated in the lower half 
of Figure 18. The computed pressure drop indicates 
streamline curvature at the pressure-surface/trailing-edge 
comer accompanied by a local flow acceleration toward 
the trailing-edge. The absence of a drop in pressure for 
the measured pressure-surface data indicates that this 
acceleration with deflection does not occur in the experi- 
mental impeller. 

The reason(s) why the numerical simulations fail to 
properly model the flow physics at the impeller exit is not 
presently known, but the authors feel confident that it is 
not due to the basic numerical scheme or approach; for 
example, grid density, Runge-Kutta time-integration, arti- 
ficial dissipation scheme (all three artificial dissipation 
schemes fail almost equally), etc. An earlier, independent 
computational study by Hathaway et al. [32] on the same 
impeller revealed similar results using a different Navier- 
Stokes code with different computational grids. It is cur- 
rently thought that some important physical effect has not 
been modeled adequately if at all, such as large-scale 
unsteadiness (vortex shedding), or that an important exper- 
imental boundary condition (related to the LSCC test rig) 
has possibly been neglected. The turbulence model may 
also be an important factor, especially as it relates to wake 
turbulence. 

In any case, the relatively large differences between 
certain computational and experimental quantities at the 
impeller exit (station 178) make it more complicated to 
use the experimental data there to confirm or discredit the 
relative accuracy of the three CFD solutions. Fortunately, 
the results presented previously seem adequate to establish 
that the ASLIP scheme is by far the most accurate, fol- 
lowed by the baseline scheme with preconditioning. With 
this in mind it can be pointed out in Figures 15 through 17 
that the more accurate the scheme, the finer the apparent 
resolution of viscous flow field features as revealed by 
spanwise undulations in the curves. These curves, how- 
ever, represent circumferentially-averaged data, so it is of 
interest to examine and compare corresponding contour 
plots, as will be done next. These comparisons will fur- 
ther solidify the conclusions about scheme relative accu- 
racy, as well as demonstrate that the impeller exit flow 


field computed using the ASLIP scheme is fairly accurate 
in many respects despite the problems delineated above. 

Contour plots of radial velocity V r /u iip at the 

impeller exit (station 178, see Figure 9) are shown in 
Figure 19a where the dashed boxes enclose the area in 
which experimental LDV data were acquired. The experi- 
mental contours are shown at the top of the figure, fol- 
lowed underneath by the ASLIP, baseline with 
preconditioning, and baseline without preconditioning 
results. Very noticeable is the increase in flow-field reso- 
lution while progressing from the non-preconditioned, to 
the preconditioned baseline, to the ASLIP solution, the lat- 
ter agreeing most closely with the LDV data. Qualita- 
tively there are some features in the ASLIP solution not 
appearing in the experimental result, for example the low 
velocity zone near the upper-left comer of the dashed box, 
but overall the agreement is remarkably good. Also, 
where there is good qualitative agreement, the quantitative 
comparison is fairly close as well. 

A similar figure is presented next, Figure 1 9b, which 
contains contour plots of relative circumferential velocity 
V e /U tip- The same general observations can also be 

made from this figure except, for reasons already dis- 
cussed, those involving direct quantitative comparisons 
between computation and experiment. However, if a con- 
stant value of about 0.09 is subtracted from the CFD con- 
tours then the quantitative agreement is also fairly good 
for the ASLIP solution. 

Finally, to complete the detailed comparisons several 
contour plots of the flow field at three locations within the 
impeller, stations 85, 126, and 165 (see Figure 9), will be 
presented. At these locations the experimental results will 
be compared only with the ASLIP results since they are 
the most accurate. 

The first comparison at station 85 is near the inlet of 
the impeller, around 15 percent of meridional distance 
behind the impeller blade leading edges. At that location 
the flow is predominately axial and just a relatively small 
fraction of the total impeller work has been done on the 
fluid. Therefore only axial velocity V X /V tip contours, 

shown in Figure 20, are compared there. In the figure, the 
ASLIP results are on the left while the experimental LDV 
results are within the dashed frame on the right. As can be 
seen, the agreement is very good qualitatively and quanti- 
tatively. 

The next comparison at station 126 is just over half- 
way through the impeller, at about 56 percent of meridi- 
onal distance from the blade leading edges. At this 
location contour plots for all three components of velocity, 

v x /u tip> v r /u tip' and v e /u ,ip • are shown in 

Figure 21, again with the ASLIP results on the left and the 
LDV results in the dashed frame on the right. The agree- 


NASA TM-1 13120 


20 



Velocity Contours @ J = 178 

Increment 0.01 

Radial Velocity, V / U fjp 


0.21 


Rotation 


LDV 

Data 


ASUP 


°-2e 


0.26 


O’ 


0.30 


0.22 






0.19 




0.32 


W 


<v 

c> 


•A. 


0.21 


0.23 


0.30 


0.27 




0.22 


Baseline 

(Preconditioning) 


0.31 \ 


0 . 2 ® 


r 

f 




0.21 


0.23 L - 


0.29 




Baseline (No 
Preconditioning) 


0.28 


$ 


0.27 


0.24 


Figure 19a — Measured and computed radial velocity contours at station 178 


ment between computation and experiment is fairly good 
for all three components, although careful examination 
does reveal some minor differences. 

The last comparison is near the exit of the impeller at 
station 165 where spanwise distributions of circumferen- 
tially-averaged velocities were presented earlier (see 
Figures 16 and 17). Contours of radial velocity V r / U ( ip 

and relative circumferential velocity V 0 V U ti are shown 

in Figure 22 in the same format used for the previous two 
figures. As can be seen, there is again close agreement 
between computation and experiment for both compo- 
nents, except within the tip-clearance vortex region 
located in the upper-right quadrant of the blade passage. 


For radial velocity the discrepancies are relatively small, 
but for circumferential velocity a very substantial differ- 
ence exists at the core of the vortex where viscous effects 
dominate. Specifically, in the computational result the rel- 
ative circumferential velocity decreases monotonically 
toward the center on the region, reaching a minimum value 
of 0.21. In the experimental result, however, the relative 
circumferential velocity first begins to decrease and then 
increases toward the center of the region which has a mea- 
sured value of about 0.46. This basic difference between 
computation and experiment was also apparent in 
Figures 16 (bottom) and 17 (bottom) at around the 80 per- 
cent span location. 
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Figure 19b — Measured and computed relative circumferential velocity contours at station 178 


Conclusions 

A preconditioning scheme has been implemented into 
a three-dimensional viscous computational fluid dynamics 
code for turbomachine blade rows. The preconditioning 
allows the code, RVC3D, originally developed for simulat- 
ing compressible flow fields, to be applied to nearly- 
incompressible, low Mach number flows. A brief descrip- 
tion was given of the compressible Navier-Stokes equa- 
tions for a rotating coordinate system, along with the 
preconditioning method employed. Details about the con- 
servative formulation of artificial dissipation were pro- 
vided, and different artificial dissipation schemes, namely 
the baseline scheme, the SLIP scheme, and the ASLIP 


scheme, were discussed and compared. A preliminary 
comparison for a laminar flat-plate boundary layer was 
made first, followed by the application of the precondi- 
tioned code to a well-documented case involving the 
NASA large low-speed centrifugal compressor for which 
detailed experimental data are available for comparison. 
The computational and experimental results, including 
performance and flow field data, were compared at the 
near-design operating point of the compressor. 

Results for three numerical schemes were investigated 
in detail, the baseline (artificial dissipation) scheme with- 
out preconditioning, the baseline scheme with precondi- 
tioning, and the ASLIP scheme with preconditioning. 
Several conclusions can be summarized from the study: 
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Figure 20— Measured and computed axial velocity con- 
tours at station 85 


• Preconditioning greatly accelerates the convergence 
rate of the numerical simulation, while also improving 
the accuracy of the steady-state solution (for the low 
Mach number range investigated). 

• The solution computed using the ASLIP scheme was 
considerably more accurate than the baseline solution 
(with preconditioning). The laminar boundary layer 
results revealed relative levels of accuracy and non- 
physical flow features (velocity profile over-shoots) 
which were representative of the more-complicated 
viscous three-dimensional solution results. A disad- 
vantage of the ASLIP scheme was poorer conver- 
gence behavior, resulting in longer solution times, but 
the improvement in accuracy is considered to be well 
worth the added computational overhead. 

• In general, good agreement was observed between the 
computed (ASLIP scheme) and experimental results. 
An area of relatively poor agreement, the impeller 
blade loading and work at the impeller exit, was iden- 
tified and addressed in some detail. It was suggested 
that the problem is possibly due to inadequate model- 
ing of some important physical effect, or to the 
neglect of an important boundary condition present in 
the experimental test rig. Deficiencies in the turbu- 
lence model were also cited as a possible cause of the 
discrepancy. 

In closing, it could be added that the LSCC impeller 
proved to be a good test case for the numerical investiga- 
tion reported herein. The medium-low Mach number lev- 


Velocity Contours @ J = 126 

Increment 0.01 


Axial Velocity , 1 / / U fj 



Radial Velocity, V/U. 

7 r tip 



Relative Circumferential Velocity , V Q 7 U fjp 



Rotation 


Figure 21 — Measured and computed axial, radial, and rela- 
tive circumferential velocity contours at 
station 126 
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Figure 22 — Measured and computed radial and relative circumferential velocity contours at station 165 


els made a comparison with the non-preconditioned 
solution practical, and the viscous-dominated flow field in 
the aft-half of the impeller exposed substantial differences 
which might not have otherwise been as apparent. 
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